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In this paper we study the process of the subsequent (runaway) fragmentation of the rotating 
isothermal Giant Molecular Cloud (GMC) complex. Our own developed Smoothed Particle Hydro- 
dynamics (SPH) gas-dynamical model successfully reproduce the observed Cloud Mass-distribution 
Function (CMF) in our Galaxy (even the differences between the inner and outer parts of our 
Galaxy). The steady state CMF is established during the collapse within a free-fall timescale 
of the GMC. We show that one of the key parameters, which defines the observed slope of the 
present day CMF, is the initial ratio of the rotational (turbulent) and gravitational energy inside 
the fragmented GMC. 



INTRODUCTION 

SPH based 3D hydrodynamical codes, starting with the series of pioneering works by Monaghan and Lattanzio 
[HI El 1131 ITTj , are always very successfully applied to the study of evolution and fragmentation in molecular 
clouds and molecular cloud complexes. These early simulations have been usually performed with a few hundred 
to a few thousand of SPH particles and with a fixed (few parsec) spatial resolution. 

Nowadays the most up to date simulations of molecular cloud evolution {e.g. |TS]) are performed using a few 
tens of thousands of SPH particles with variable smoothing lengths. These simulations also include the details 
of cooling and heating in the complex gas mixtures of H, H2, CO and HII species. 

Our present high resolution (64,000 SPH particle) simulations, with highly flexible and adaptive smoothing 
lengths, study the runaway collapse and the subsequent isothermal fragmentation of the isolated GMC complex 
with different rotational (turbulent) energy parameters of the clouds. The resulting CMF is compared with the 
recent observational distributions (c?N/dM ~ M r , where the slope of the power law 7 is in the range 1.4 to 
1.8) of the molecular cloud complexes derived from the different CO data for the different parts of our Galaxy 

[ni [23 m in mi Em 



METHOD 

Continuous hydrodynamic fields in SPH are described by the interpolation functions constructed from the 
known values of these functions at randomly positioned N "smooth" particles with individual masses |18| . 
To achieve the same level of accuracy for all points in the fluid it is necessary to use a spatially variable smoothing 
length. In this case each particle has an individual value of the smoothing length - hi. 

A more detailed and complete description of the basic numerical equations of SPH can be found in many of 
our previous publications (e.g. [Tl 121 131 El I24j and the references herein). Therefore we just briefly repeat the 
skeleton SPH equations of the code here. The density at the position of the particle i can be defined as: 

N 

Pi = ^ "' j • Wij, 

3 = 1 
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Figure 1. Time evolution of the thermal, kinetic, potential and total energy for the adiabatic collapse of an 
initially "cold" gas sphere. The different lines corresponds to the different gas particle numbers. We also present 
the energy error results of the GADGET code with "standard" parameters for 32,000 particles. 



The equations of motion for a particle i: 



~dt 



r/v IP P \ 



where Pi is the pressure, $i is the self gravitational potential and Ely is an artificial viscosity term. 
The internal energy equation has the form: 



du t 1 A (P t Pj - \ 



Pi 



Here Ui is the specific internal energy of the particle i. The term (r^ — hi)/ pi accounts for non adiabatic 
processes not associated with the artificial viscosity. We present the radiative cooling in the form proposed by 
[T5| (see case "B") using the MAPPINGS III software 

A = h(p, u, Z, ...) ~ h*(T, [Fe/H] ) • nj, n, = Pi /(p ■ m p ), 

where m is the hydrogen number density, Ti the temperature and p the molecular weight. 
The equation of state must be added to close the system: 



Pi = il - 1) Pi ■ Ui, 

where 7 is the adiabatic index. 

In SPH one of the basic tasks is to find the nearest neighbors of each SPH particle, i.e. to construct an 
interaction list for each particles. Basically we need to find all particles with | |< 2max(/i.;, hj) in order to 
estimate the density and also calculate the hydrodynamical forces. 

In our code we keep the number of neighbors exactly constant by defining 2hi to be the distance to the Nb 
- nearest particle. The value of Nb is chosen such that a certain fraction of the total number of "gas" particles 
N affects the local flow characteristics. From these we need to SELECT the closest Nb particles. Fast algorithms 
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Figure 2. The time evolution of the total number of clusters N c and the total mass fraction inside these clusters 
<p during the simulations. Starting from the ~ 5 Myrs more than 80 % of the total mass is already concentrated 
inside the fragments. At around 6 Myrs already almost 95 % of the total mass is inside the clusters. 



for doing this exist [201 • For computational reasons, if the defined hi becomes smaller than the selected minimal 
smoothing length h m im we set the value hi = h m i n . 

To calculate the self gravitational potential and self gravitational force — Vi$i we use the Mitaka Un- 
derground Vineyard (MUV) GRAPE6 computer system at the National Astronomical Observatory of Japan 
http : //www . cc . nao . ac . jp/muv/ . For a more detailed description of the GRAPE6 board and fot links to pub- 
lication about the GRAPE6, we refer the reader to the official homepage of Jun Makino at Tokyo University 
http : //grape . astron. s .u-tokyo . ac . jp/~makino/grape6 .html . 

For the time integration of the system of hydrodynamical equations we use the second order Runge-Kutta- 
Fchlberg scheme. The time step Ati for each particle depends on the particle's acceleration and velocity Vj, 
as well as on the sound speed Cj and the heating vs. cooling balance: 



At = C n 



2hj hj tn V4 

I a,; I ' I V,; I ' Ci ' Ui 



where C n is the Courant's number = 0.1. For computational reasons we fix the minimal integration time step 
At m ^ n . 

The main aim of our current work is a detailed study of the isothermal fragmentation processes inside 
the collapsing "cold" molecular cloud complex. For the purpose of finding the fragments and its physical 
parameters (mass and size) we use our own cluster finding algorithms. In our algorithms we modify the well 
known and "standard" friend-of- friend (FOF) method (TJ. Instead of just using the particle positions in the 
process of "constructing" or finding the clusters (fragments) we also use the information about the density 
distribution inside each potential cluster. On the basis of the density distribution analysis we can finally select 
in the more accurate way the members of our fragments (clusters). In this sense our method is more close 
to the so called SKID method, which is well described at the homepage of the Washington University "N- 
body Shop" http://www-hpcc.astro.washington.edu/tools/skid.html . Here the reader can find a more 
detailed description of this density base method which is specially designed to find the gravitationally bound 
groups of particles in the N-body like simulations. 

One of the features of our cluster finding routine is in the setting of the minimum limit of gaseous particles 
to 5, in order for a fragment to form. In other words we don't count as "real" a cluster where the number of 
particles is less then 5. 
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c. "Fast" initial rotation uo = 0.3, Time 
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b. "Slow" initial rotation w = 0.1, Time = 7 Myr. 
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Figure 3. Four snapshots of the integrated cluster distribution function (ICMF) 



CODE TESTING 

The selfgravitating collapse of an initially isothermal "cold" gas sphere has been a common test problem for 
different SPH codes [7| ES] |3J [Z7\ 123 • Following these authors, for the testing of our code, we calculate the 
adiabatic evolution of the spherically symmetric gas cloud of total mass M and radius R. For the initial internal 
energy per unit mass we set the value: u = 0.05 ^jr^- The initial density profile of the cloud calculates as: 

( \ M 1 

2 7r Hr r 

We distribute randomly the gas particles inside the set of spherical shells in a manner that reproduces the 
initial density profile. At the start of the simulation the gas particles are at rest. For the presentation of the 
results we use a system of units where G = M = R = 1. 

In Fig. n we show the time evolution of the different types of energy and the relative total energy error 
during the calculation. For comparison of our test results we also plot the energy error results from the serial 
variant of the GADGET public access SPH-TREE code [231 with "standard" parameters for the 32,000 particles 
http : // www . mpa-gar ching . mpg . de/gadget / . 
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Figure 4. The ICMF slope time evolution for our models with different rotation parameters. 

During the central bounce around i w 1.1 most of the kinetic energy is converted into heat, and a strong 
shock wave travels outward. For all of these runs the number of neighbors was set Nb — 50 and the 
gravitational softening was set e — 0.01. For the integration of the system of equations we use the second 
order Runge-Kutta-Fehlberg scheme with a fixed time step At = 10 -4 . 

The results presented in Fig.^agree very well with those of j2U and |23|. The maximum relative total energy 
error is around 0.05 % even for moderate (8,000) particle numbers. The largest adiabatic test calculation (with 
64,000 gas particles up to t rj 2.2) on an Intel Pentium 4 (3.4 GHz) host machine with a GRAPE6 board took « 
3.67 days of total CPU time. 

INITIAL CONDITIONS 

As an initial condition for our molecular cloud fragmentation study we use a model in which the parameters 
are comparable with the largest GMC complexes in our Galaxy ^1 EH QUI ■ For the mass of the system we set 
M c i OU( 2 = 10 7 Mq. For the radius of the cloud we set R c i ou d = 100 pc. For an initial density distribution we 
use the previous formula where p{r) ~ K For the purpose of checking the possible "resolution" effects we carry 
out two sets of runs with 32,000 and 64,000 gas particles (with the corresponding indexes "low" and "high"). 
The total gravitational energy of the system in such a case can be easily calculated using the simple formula: 

Eqra — — g G M cloud /R c i oud . 

For the initial temperature we set the value which produced the overall ratio of the thermal energy to 
the gravitational energy of the system at the fixed level a = E^ he /|EqraI = 0.075. For the previous fixed 
mass and radius of the system this condition produced an initial temperature of the cloud T c ioud ~ 2200 K. 
The corresponding sound speed was c « 3.8 km/sec. This is consistent with the typical measured "kinetic" 
temperatures for such GMC complexes UJ. 

With these parameters we have an initial central concentration of no w 10 3 cm -3 , and a free-fall time in the 
cloud center of Tfj « 1 Myr. The central Jeans radius is Rj « 10 pc with the corresponding Jeans mass of Mj 
« 10 5 Mq. Initially we give the whole system a rigid rotational velocity distribution with an angular velocity 
value £l c ioud which we set to the unity of Qq = Vq/ RcZoud where: 
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V = M cloud /R 

cloud- 

Using our parameters this velocity is equal to Vo ~ 21 km/sec. 

The main rotational parameters (lo and (3) for our two sets of models are listed in Tabled 



Table 1. The list of initial "rotational" parameters in our models 
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= 0.05 
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« 6.2 10" 4 


2.5 10" 3 


5.6 10" 3 


1.0 10" 2 


1.5 10" 2 


2.2 10" 2 



As we can see from Table H the initial ratio of the rotational (or kinetic) energy of the motion the the 
fragments (/3), even in a last models with "high" rotational parameter, don't exceed more than a few % from 
the gravitational bounding energy of the cloud. This is even less for the initial ratio of the systems thermal 
energy to the gravitational energy (a = 7.5 %). In all models, usually after the first few Myr of evolution, these 
situations have changed. The ratio (3 is rising to the approximate value of 0.5 or even more. Its mean is when 
the cloud starts the process of intensive isothermal fragmentation and the whole system of fragments becomes 
almost fully "rotationally" supported. 

RESULTS 

In Fig. [5] we show the time evolution of the total number of clusters N c and the total mass fraction inside 
these clusters <fi during the simulations. Starting from ~ 5 Myr more than 80 % of the total mass is already 
concentrated inside the fragments. At around 6 Myr already almost 95 % of the total mass is inside the clusters. 
Fig. Et shows the results for the "slow" rotating model with lo = 0.1. In Fig. Eh we show the evolution of the 
"fast" rotating model with lo = 0.3. 

In Fig. 01 we show four different snapshots of the integrated cluster distribution function (ICMF) for two 
selected moments of time with two different rotational parameters lo. For practical numerical reasons we use 
the ICMF instead of the CMF (which is sometimes in the literature also called the Differential Cloud Mass 
Function DCMF). Here is a simple definition of the ICMF: 

ICMF= / dN/dM. 
Jo 

Basically it shows how many clouds we have from zero mass to any fixed mass (M). Because the CMF is 
usually approximated with the power law: M -7 in this case the ICMF we can be simply derived as ~ M _7+1 . 

The reason for using the ICMF instead of the CMF (DCMF) is that the averaging and slope definition is 
mathematically better due to the integrated CMF (which is a monotone function) and because the histograms 
in this case don't have any "holes" . Of course when we compare our results with the observed (differential) 
CMF slope we need to subtract one from the ICMF slope to get the corresponding CMF slope. 

In Fig. iniwe can see that in most cases the ICMF slope lies between -0.5 and -0.7 (the corresponding CMF 
is -1.5 and -1.7). The models with slow rotation always have a significantly lower value of the slope. 

The ICMF slope time evolution for the set of our models with different rotation parameters are presented in 
Fig- El The models with initial "slow" rotational parameters give the ICMF average slope a level of -0.8 (which 
corresponds to the CMF slope -1.8). The "fast" rotating models give the ICMF slope a level of -0.4 (CMF slope 
-1.4). 

On average all models show very close values of the CMF slopes in comparison with the observed values of 
7 in the different parts of our Galaxy. 

The slow rotation models systematically show the slope more close to the observed values in the outer part 
of our Galaxy ^1^3 |S] (7 « -1.8±0.03). In contrast, the fast model CMF slopes is more consistent with the 
observations from the central part of our Galaxy 7 ~ -1.5 ±0.1. 

Of course our simulations are time and also resolution limited, but even in this case we can derive a state- 
ment about the two significantly different types of "population" in the molecular cloud distributions. The key 
parameter which produces the different CMF slopes is the initial rotational parameter of the forming (and 
subsequently fragmenting) GMC. 
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CONCLUSIONS 



In this paper we present a study of the subsequent (runaway) fragmentation of the rotating isothermal GMC 
complex. Our own developed GRAPE based Smoothed Particle Hydrodynamics (SPH) gas-dynamical model 
successfully reproduced the observed Cloud Mass-distribution Function (CMF) in our Galaxy. The steady state 
CMF is quickly established during the collapse approximately on a scale of a few free-fall time in the central 
parts of the modeled GMC. 

One of the key points in our model is that using our results we can naturally explain the source of possible 
differences between the observed slope on molecular clouds mass distribution function in the Galactic center 
and the outer regions of our Galaxy. 

The basic idea, is what if the GMC formed as a result of the galactic disk instability on the scale of the 
disk height (~100 pc). In such a case the initial angular momentum of the forming GMC can be defined by 
the Coriolis force during the formation inside the differentially rotating disk. Therefore the central GMC has a 
bigger (3 and the external GMC has a smaller rotational parameter. 

According to our models this produces the different slopes of the resulting CMF during the runaway frag- 
mentation process inside the system. The observed CMF gives to the central parts of Galaxy a slope well 
approximated with the value 7 « -1.5±0.1 |221 1211 [TB] and for the outer parts of the Galaxy the approximate 
value 7 « -1.8±0.03 H2H333E1. 

Our results for the "slow" and "fast" rotating models give us exactly the same slopes with very good 
agreement with the recent observations. The "slow" models corresponds to the initially more slowly rotating 
GMC in the outer parts of the Galaxy. The "fast" rotating models corresponds to the GMC in the central part 
of the Galaxy. The central GMC can initially get more angular momentum from the differential rotation of the 
galactic disk during the process of GMC formation itself. 

Our numerical investigation clearly shows that one of the key parameters, which determines the observed 
slope of the present day molecular CMF in different parts of our Galaxy, is the initial ratio of the rotational 
(turbulent) and gravitational energy inside the forming GMC. 
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